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Abstract. We propose a high order numerical decomposition of exponentials of 
hermitean operators in terms of a product of exponentials of simple terms, following 



^•J^ ■ an idea which has been pioneered by M. Suzuki, however implementing it for complex 

C^) ' coefficients. We outline a convenient fourth order formula which can be written 

' compactly for arbitrary number of noncommuting terms in the Hamiltonian and which 

| is superiour to the optimal formula with real coefficients, both in complexity and 

*^ accuracy. We show asymptotic stability of our method for sufficiently small time step 

Oh, and demonstrate its efficiency and accuracy in different numerical models. 

cr. 



1. Introduction 

While exponentials of operators are very common in every field of quantum physics, but 
also in classical physics, their evaluation is nevertheless numerically a very demanding 
operation. For example, in quantum physics, this task usually emerges when one wants 
to compute a time-evolution, either in real time, for example when computing dynamical 
correlations, or in imaginary time, when computing thermodynamic averages like in 
quantum Monte Carlo simulations. A similar decomposition of classical time evolution, 
which can also be interpreted in terms of unitary operators, is known as symplectic 
integration. 

For an operator which can be written as a sum of several parts of which exponential 
operators are exactly determinable, the well known Suzuki- Trotter [TJ 121 El IH EH EH E] 
decomposition scheme can be used. The operator e lz ^ Aj is approximated by a product 
of operators e lzpfc j Aj with real coefficients such that the desired order of accuracy is 
achieved. We will show in the present paper that following the same principles but not 
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restricting to real coefficients the same order can be achieved using a smaller number 
of factors. Furthermore, the order of such decomposition can be trivially increased by 
one by composing it with an equivalent decomposition with a complex conjugate set 
of coefficients. We will outline a particular third order scheme, and further improved 
to fourth order, which is potentially very useful for practical calculations. We show 
explicitly that, even though we lose unitarity of decomposition (in real-time case), the 
method is asymptotically stable for sufficiently small time steps since all the eigenvalues 
of the decomposition remain on the complex unit circle. Even more generally, we show 
that one gains an extra order in accuracy and asymptotic stability (independent of the 
size of the time step) by renormalizing the state vector after each time step. 

We demostrate the accuracy and efficiency of the method by three explicit examples: 

(i) in case of 2 x 2 matrices the decomposition and its stability can be treated analytically, 

(ii) for exponentials of Gaussian random Hermitean matrices we find that the stability 
threshold (the maximal time-step for which the method is asymptotically stable) drops 
with the inverse power of the dimension of the matrix, and (iii) for a generic (non- 
integrable) interacting spin 1/2 chain (in one-dimension) we find, surprisingly, that the 
stability threshold is independent of the number of spins. 

2. Complex Split-Step Decomposition 

Our main objective is to approximate the exponential operator Uq = e lz ^ A+B \ for general 
bounded operators A and B, and a complex parameter z, as a product U of exponential 
operators 



The equations determining the coefficients {pj} that solve the equation above are 
obtained by expanding the exponential operators into power series and equating lowest 
order terms to zero. It is known that there is no third order (0(z 4 )) solution of the five- 
term ansatz (JIJ with real coeffients pj. The simplest third order decomposition involves 
six terms [Hj. However, allowing the coefficients pj to be complex, there exist two very 




(1) 
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simple and symmetric solutions, namely X 

_ 1 V3. _ 1 y/3. 1 

Pi = P5 = 4 + j2 i, P2 = p 4 = g + "g" 1 . ^3 = 2 ( 2 ) 

and the complex conjugate set {pj}. 

Let us denote exact exponential as Uq(z) = exp(iz(A + B)) and third order complex 
decompositions (C3), given by RHS of (JTJ) with coefficients (J2J), namely {pj}, and [p]}, 
as £/(;z) and £/(<z)> respectively. Using some further analysis (which has been performed 
by means of Mathematica software) we can show that the next-order-term changes sign 
when one switches between the two solutions, namely: 

U(z) = U (z) + K A z 4 + 0(z 5 ) and U(z) = U {z) - K 4 z* + 0{z 5 ), (3) 

where 

K A = ^ ((AAAB - BAAA) - 3(AABA - ABAA) - 
- 3(AABB - BBAA) + Q(ABAB - BAB A) + 
+ 2(ABBB - BBBA) + Q(BABB - BBAB)) (4) 

is a Hermitean operator provided that both A and B are Hermitean. 

Superposition of the two decompositions cancels the term and is therefore for one 
order higher, namely of fourth order. However, the same, fourth, order can be achieved 
by alternating both decompositions (as illustrated in fig. ^) 

U(z)U(z) = Ul + (U K 4 - K 4 U )z 4 + 0(z 5 ) = Ul + 0(z 5 ), (5) 

since Uq(z) = l+0(z). Since in usual numerical simulations of exponential operators, for 
example in quantum time-evolutions, time dependent renormalization group methods, 
or quantum Monte-Carlo simulations, one needs to make many time-steps anyway, the 
alternation between U(z) and U(z) does not represent any practical drawback. 

However, we note that with pi being complex numbers the decomposition U(z) is 
no longer strictly unitary (in the usual case where the operators A and B are Hermitean 
and the time step z is real) and the time evolved state (on which U operates) might 
explode in norm after a while. In order to strictly preserve the norm, the state (vector) 
may be renormalized at every time step. One might be afraid that this renormalization 
X It was quoted in Ref. 6 that this solution had already been proposed by A.D.Bandrauk, however it 
was claimed in Ref. 8 that the complex coefficient decomposition is unstable and cannot be practically 
used for splitting the unitary exponentials, which we show is not precise. 
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Figure 1. Schematic illustration of complex valued split step decomposition. 
Coefficients pj can be considered as shifts in complex time plane, which always 
move along the real axis. Two sets of complex coefficients {p{\ give a third order 
decomposition C(z 4 ); their superposition is for an order higher. 

would degrade the accuracy of the method. However, due to the fact K\ = K4 this is 
not the case, in fact renormalization increases the accuracy to fourth order 

(Uj{z)U{z)) _ l + (K 4 )z* + Q(z 5 ) _ 1|g( , 5) (6) 

yJ{W(z)U(z)) yJl + (Ki + Kl)z* + 0(zS) 

By (.) := (^l-l^) we denote the expectation value in some intial state vector In 
conclusion, the decomposition with one single set of complex coefficients pi (C3) is 
already of the fourth order (0(z 5 )) if every time step is followed by renormalization of 
the state (fig. EJ). As in any application the computational complexity of performing 
the sequence of exponential operators on a state vector U(z)\tp) is dominating the 
normalization of the state, this does not represent any drawback of the method. Still, 
as we will show later, the method is asymptotically stable, for sufficiently small z even 
without the renormalization. Figure |2] shows real numerical errors, in a model in which 
A and B are chosen as Gaussian random Hermitean matrices, after performing two 
time steps with various decompositions described above (using one (C3) or both sets 
of complex coefficients (C4), and with or without renormalization of the state) and 
compare it with the optimal third order decomposition with real coefficients (R3). 

We can easily generalize our approach to approximate exponentials of three or more 
noncommuting bound operators. For example, for three operators, one has nine terms 
following a sequence ABC B ABC B A which is obtained from ABABA ((TJ by replacing 
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Figure 2. An error after two time steps for the third-order real decomposition 
(R3), the third-order complex decomposition (C3), and the fourth-order complex 
decomposition (C4); the label V denotes renormalization after each time step. As 
for numerical model we choose A and B to be GUE matrices of dimension A = 200 
and average the results over 1000 realizations. 



each inner operator B by BCB (and dividing the coefficient in front of B by two) 



Az(A+B+C) 



Azp 1 A\zp 1 B \zp 2 C \zp\B izp^A izp4,B \zp 5 C \zpiB izp 5 A 



and using the same set of coefficients (|2|L or its complex conjugate. Generally, a formula 
for a sum of n operators involves An — 3 terms 

exp (iz(At + . . . A n )) = 



izpiAi Azp\A 2 izpiAn-i AzpiA n Azp\A n -\ a \zp\A 2 v 

e e • • • e e e • • • e a 

^zp-iAx e izp 5 A 2 _ _ _ Q izp 5 A n - 1 ^izp4,A rl ^izp 5 A n - 1 _ _ _ ^-zp 5 A 2 ^\zp 5 A 1 



(8) 



It is interesting to note that the general optimal third order solution with real coefficients 
(R3) uses just one term more for the case n = 2, namely six, whereas for general n case 
it needs 5n — 4 terms, which is n — 1 terms more than the complex solution above (JHJ). 

As we have mentioned before, without the renormalization complexness of the 
coefficients may cause the exponential instability of the method. However, it turns 
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out that the decomposition is absolutely stable for small enough steps z. The reason 
for such an interesting behaviour is that the eigenvalues of the operator U(z) lie all 
on complex unit circle for sufficiently small z, and this property grants the asymptotic 
stability even if U(z) is not exactly unitary There is typically a threshold, i.e. a critical 
value of z max such that at z = z max two eigenvalues of U(z) collide and leave the unit 
circle and then the method ceases to be asymptotically stable. Such a behaviour can be 
explicitly proven for operators chosen from the space of 2 x 2 matrices (see the following 
section) and is conjectured in general. 

3. Examples 

First, let us consider a numerical example of calculating the exponential of H = A + B 
where A and B are Gaussian random Hermitean matrices chosen at random from the 
Gaussian Unitary Ensemble Figure EK shows that the maximal size of eigenvalue of 
U(z) is exactly equal to one until some point described by the threshold step size z max - 
Numerical results suggest the following dependence of the threshold on the Hilbert space 
dimension N, z max oc 1/N a , with a ~ 0.5, which we believe is the worst case scenario 
for generic systems. 

As a second example, we consider a non-trivial physical model where the matrices of op- 
erators A and B are very sparse and thus far from the full random matrix model, namely 
we consider time evolution in the quantum Ising spin 1/2 chain in a tilted homogeneous 
magnetic field (e.g. recently considered in the context of heat transport JUj) described 
by the hamiltonian H = £n=i + g x o% + 9zO" z n }. Here, cx^' z , n = 1 . . . N s , 

represent a set of independent Pauli matrices. In figure Efc> we show a very interesting 
result for this model (in particular, for the parameter values J — 1, g x = 0.4, g z = 0.8 
which lie in the so-called "quantum chaotic" regime JU]), namely that the threshold 
step size z max is asymptotically independent of the size iV = 2^ of the system. We 
conjecture that this is in general true for numerical simulations of finite (spin) quantum 
systems with local interaction, and for such our method of simulation of time-evolution 
should be very roboust. 

As for the last example, we make analytical consideration of the simplest case where our 
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Figure 3. Maximal size of the eigenvalue of the approximate evolution operator U (z). 
The upper plot (a) shows the case of GUE matrices while the lower plot (b) shows 
the case of Ising spin chain in tilted magnetic field (see text). Different curves refer 
to systems of different sizes (b), or different matrix dimensions (a). The insets show 
critical threshold z max as a function of the system size/matrix dimension. 
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operators can be represented by 2 x 2 matrices. In order to understand the transition 
in the stability (collision of eigenvalues of U(z) on the unit circle) one can generally 
parametrize the operators A and B by Pauli operators a\ j = 1, 2, 3, 

3 3 

A = a l + and B = b l + ^ b j( r j . (9) 

The coefficients {cl,-}, and {bj} are all real since matrices A and B are Hermitean, and 
furthermore matrices A and B can always be chosen traceless by setting ao = bo = 
without losing generality. It is obvious that, since det U = e izTrH , where H = A + B, 
that decomposition for two 2x2 matrices can also be expressed in terms of Pauli 
matrices and some coefficients {gj}- Using the ansatz (0) we write 

Of course, gj are no longer real in general. Eigenvalues of the operator U (z) = e 12 ^ 9](j3 



are e V 3 3 which gives the condition for the asymptotic stability: namely the number 
7 2 = J2j Q 2 j should be real and positive, j 2 G R + . In order to simplify the notation, let 



us take 7 = +yJ2j9j> an d similarly write a = yJ2j a p(3 = yY^jbj, and introduce 
normalized coefficients ji = gi/j,aj = aj/a,/3j = bj/j3. The condition for asymptotic 
stability now simply reads 7 G R. Using straightforward calculation 7 can be expressed 
as 7 = - arccos (iTre^ 9 ^) and is, interestingly, only a function of the magnitudes a, 
(3 and z-projections a 3 and /? 3 : 

7(2) = - arccosQ(^), where 

z 

Q(z) = ^(l - a§ + (1 + 3a 2 ) cos(az)) ((1 + /3 3 2 ) cos(/3z) + 



+ (1 - cosh(-^)j - 2a 3 (3 + a 2 3 )f3 3 sin(az) sin(f3z) + 
V3 



+ 2(1 - a 2 ) cosh(^) ((1 + (3 2 3 ) cos(^) cos(^) + 

+ (l-/5 3 2 )cos(^)cosh(^|)-2a 3 /5 3 sin(^)sin(^))^ (11) 

Now the stability condition reduces to |Q(^)| < 1. For small steps z the expression Q(z) 
in (fTTj) can be written as a power series in z 

Q(z) = 1 - Ua 2 + (3 2 + 2a 2 + 26 2 + 6a 3 6 3 )z 2 + 0{z^). (12) 


It can easily be proven diagonalizing the quadratic form that the z 2 term is always 
nonpositive, hence the decomposition scheme indeed is always (for any a,j, bj) stable, for 
small steps z. 
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Figure 4. Illustration of the stability threshold for 2x2 case. Since matrices are 
traceless, collision of eigenvalues of U (z) takes place on the real axis. In the figure 
we plot Re7 (dashed) and Inry (full), as a function of z for the case a = j3 = 1 and 
a 3 = h = 0.1. 

Figure H] illustrates how eigenvalues for small steps z always lie on the unit circle 
in the complex plane. When the step z is being increased, the eigenvalues are travelling 
along the unit circle, one in clockwise and the other in the counter-clockwise direction. 
At some point, namely at z = z max , a collision occurs and a pair of eigenvalues bounce off 
the unit circle - then 7 becomes complex. However, because of the restriction |det U\ — 1 
their product remains on the unit circle. Our 2x2 matrices A and B are assumed to 
be traceless therefore collisions always occur on the real axis and eigenvalues are both 
real during the bounce. 

4. Conclusion 

We have proposed a simplex explicit complex-coefficient split-step decomposition of 
an operator exponential, based on Suzuki's scheme, for a sum of arbitrary number of 
operators. As compared to an optimal scheme with real coefficients our scheme requires 
less terms for the same order, furthermore we can gain an extra order at no additional 
expense. Despite having complex coefficients the decomposition is always stable for 
sufficiently small step size, and can be stablilized by additional renormalization of the 
state vector. 
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We suggest that our method may be used in conjunction with other methods for 
efficient time evolution of complex quantum systems (one application has already been 
done in Ref.|10j). or interacting many body quantum systems, like for example with 
time-dependent DMRG methods ^lj ^] where efficient and accurate estimation of 
operator exponentials for short time steps is one of the cruicial black-box operations. 
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